Load packages

# Check you have them and load them
list.of.packages <- c("kableExtra", "tidyr", "ggplot2", "gridExtra", "dplyr", "Hmsc", "jtools", "lubridate", "corrplot", "MuMIn","stringr","sf","raster","leaflet", "tidyverse","htmlwidgets","webshot", "purrr", "magick","forcats","multcomp", "reshape2", "lme4", "tidyr", "stats","RColorBrewer")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages, require, character.only = TRUE)
## Loading required package: kableExtra
## Loading required package: tidyr
## Loading required package: ggplot2
## Loading required package: gridExtra
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: Hmsc
## Loading required package: coda
## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0
## Loading required package: jtools
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## Loading required package: corrplot
## corrplot 0.92 loaded
## Loading required package: MuMIn
## Loading required package: stringr
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.3.2
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
## Loading required package: raster
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: leaflet
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.4
## ✔ purrr   1.0.1     ✔ tibble  3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine()    masks gridExtra::combine()
## ✖ raster::extract()   masks tidyr::extract()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ raster::select()    masks dplyr::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: htmlwidgets
## 
## Loading required package: webshot
## 
## Loading required package: magick
## 
## Linking to ImageMagick 6.9.12.98
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
## 
## Loading required package: multcomp
## Warning: package 'multcomp' was built under R version 4.3.3
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.3.3
## 
## Attaching package: 'mvtnorm'
## 
## The following object is masked from 'package:jtools':
## 
##     standardize
## 
## Loading required package: survival
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 4.3.3
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following objects are masked from 'package:raster':
## 
##     area, select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## Attaching package: 'TH.data'
## 
## The following object is masked from 'package:MASS':
## 
##     geyser
## 
## Loading required package: reshape2
## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths
## 
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 4.3.2
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.3.2
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## Attaching package: 'lme4'
## 
## The following object is masked from 'package:raster':
## 
##     getData
## 
## Loading required package: RColorBrewer
## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
## 
## [[9]]
## [1] TRUE
## 
## [[10]]
## [1] TRUE
## 
## [[11]]
## [1] TRUE
## 
## [[12]]
## [1] TRUE
## 
## [[13]]
## [1] TRUE
## 
## [[14]]
## [1] TRUE
## 
## [[15]]
## [1] TRUE
## 
## [[16]]
## [1] TRUE
## 
## [[17]]
## [1] TRUE
## 
## [[18]]
## [1] TRUE
## 
## [[19]]
## [1] TRUE
## 
## [[20]]
## [1] TRUE
## 
## [[21]]
## [1] TRUE
## 
## [[22]]
## [1] TRUE
## 
## [[23]]
## [1] TRUE
## 
## [[24]]
## [1] TRUE
## 
## [[25]]
## [1] TRUE
## 
## [[26]]
## [1] TRUE

Prepare Data

locs <- read.csv("Processed_Data/TDN_camera_locations_and_covariates_cleaned.csv", header=T)

# Convert to categorical factors
locs <- locs %>% 
            mutate_if(is.character,as.factor)

# You should also standardize your covariates - it helps models coverage and facillitates comparison of effects sizes

#NEED TO RUN THIS SO LAT/LONG DONT GET STANDARDIZED
#non response variable species should be scaled but that can manually be done in the GLMM formulas
library(MuMIn)
z_locs <- stdize(locs, omit.cols = c("latitude","longitude","latitude_site","longitude_site","spatial_scale"))

Delineate Caribou Seasons

Rather than defining seasons in a community aspect using snowcover (which would be more suited to a multispecies analysis) I want to divide data by biologically relevant “caribou seasons”. The problem here is that these seasons vary by date from year to year and from herd to herd based on GPS collar data. I have 2021-2022 dates for the Ahiak herd from GN reports but I don’t have similar recent seasonal date ranges for the Bathurst/Beverly herds (date ranges used in many reports for the Bathurst herd are from a 2011 study). Do have telemetry data for the Bev/Bath caribou but it’s clipped to the spatial extent of TDN.

Determine first obs, last obs, and number of days observed in TDN for each individual collared caribou using 2021/2022 collar data. Extracted attribute table as csv from shapefile using arcgis

#read in collar data
TDN_Bev_Bath_GPS<-read.csv("Raw_Data/TDN_Bev_Bath_Telemetry.csv")

# Convert LoctnDt to Date format
TDN_Bev_Bath_GPS$LoctnDt <- as.Date(TDN_Bev_Bath_GPS$LoctnDt)

Summarize the collar data

# Group by AnimalId, Sex, and Herd, then summarize to find first date, last date, count of observations, and number of unique dates
caribou_indv_summary <- TDN_Bev_Bath_GPS %>%
  group_by(AnimlId, SEX, Herd) %>%
  summarize(
    FirstDate = min(LoctnDt),
    LastDate = max(LoctnDt),
    NumObservations = n(),
    NumUniqueDates = n_distinct(LoctnDt)
  )
## `summarise()` has grouped output by 'AnimlId', 'SEX'. You can override using
## the `.groups` argument.
# Group by Herd, Sex, and then summarize to find first date, last date, count of observations, number of individuals, and number of unique dates
herd_summary <- TDN_Bev_Bath_GPS %>%
  group_by(Herd, SEX) %>%
  summarize(
    FirstDate = min(LoctnDt),
    LastDate = max(LoctnDt),
    NumObservations = n(),
    NumIndividuals = n_distinct(AnimlId),
    NumUniqueDates = n_distinct(LoctnDt)
  )
## `summarise()` has grouped output by 'Herd'. You can override using the
## `.groups` argument.
print(herd_summary)
## # A tibble: 5 × 7
## # Groups:   Herd [3]
##   Herd       SEX    FirstDate  LastDate   NumObservations NumIndividuals
##   <chr>      <chr>  <date>     <date>               <int>          <int>
## 1 Bathurst   Female 2021-10-27 2021-11-15              45              2
## 2 Bathurst   Male   2021-10-29 2022-05-06             186              3
## 3 Beverly    Female 2021-10-27 2022-01-17             707             16
## 4 Beverly    Male   2021-10-29 2022-08-30            1049             12
## 5 Unassigned Female 2021-11-18 2021-11-21               7              1
## # ℹ 1 more variable: NumUniqueDates <int>
# If you want separate summaries for each Herd irrespective of sex
herd_summary_no_sex <- TDN_Bev_Bath_GPS %>%
  group_by(Herd) %>%
  summarize(
    FirstDate = min(LoctnDt),
    LastDate = max(LoctnDt),
    NumObservations = n(),
    NumIndividuals = n_distinct(AnimlId),
    NumUniqueDates = n_distinct(LoctnDt)
  )

Herd SEX FirstDate LastDate NumObservations NumIndividuals NumUniqueDates 1 Bathurst Female 2021-10-27 2021-11-15 45 2 18 2 Bathurst Male 2021-10-29 2022-05-06 186 3 41 3 Beverly Female 2021-10-27 2022-01-17 707 16 83 4 Beverly Male 2021-10-29 2022-08-30 1049 12 215 5 Unassigned Female 2021-11-18 2021-11-21 7 1 3

Looks like the cows are only using TDN during Rut/Breeding, Fall Migration - post-breeding, and the start of Winter according to seasonal cutoff dates other Bathurst papers have used.

Ahiak have slightly different seasonal cutoff dates but collared cows roughly use TDN in 2021/2022 in Late Summer, Fall migration - pre-breeding, Rut/breeding, Fall migration - post-breeding, Winter, no spring migration, then back for Summer and Late Summer.

If we used these dates, how do we account for the different herds (Ahiak winters here but they’re expected to winter on the tundra as opposed to Bev/Bath) and the seasons where we’re less likely to have female caribou (i.e. spring)

Multiseason, DAILY

The telemetry data was informative but still found it difficult to set seasonal cutoff dates based on that alone. Plotting daily detection rate for caribou might be a better way to find natural seasonal breaks and is the most detailed data we have on hand.

Let’s see how this breaks down into 1 day events to help determine seasonal breaks

#Load biweekly detection data
daily_obs<-read.csv("Processed_Data/wildRtrax/TDN_cam_summaries_day.csv")

Join this dataframe with the standardized location data:

mod_dat_daily <- left_join(daily_obs, z_locs)
## Joining with `by = join_by(location)`

And let’s do another raw data check:

boxplot(mod_dat_daily$Caribou~mod_dat_daily$dom_habitat_lcc_hexagon,
        las=2,
        xlab="lcc landcover",
        ylab="Habitat use")

##Add snow based seasons

#format day as Date
mod_dat_daily$day<-as.Date(mod_dat_daily$day)

#create column season
mod_dat_daily <- mod_dat_daily %>%
  mutate(snow_season = ifelse(month(day) %in% c(5, 6, 7, 8, 9, 10), "summer", "winter"))

# Get the column names
column_names <- names(mod_dat_daily)

# Reorder the columns
new_column_order <- c(column_names[1:4], "snow_season", column_names[5:length(column_names)])

# Create the reordered data frame
mod_dat_daily <- mod_dat_daily[, new_column_order]

# make it a factor factors
mod_dat_daily <- mod_dat_daily %>% 
            mutate_if(is.character,as.factor)

let’s plot caribou detections summed across locations per day coloured by snow based seasons

# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
  group_by(year = year(day), day, snow_season) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.
# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Locations by Day",
       x = "Day",
       y = "Total Caribou Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))+
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")

# Plot the data for capture rate
daily_sum$total_daily_cr <- (daily_sum$total_caribou / daily_sum$total_effort)

ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Daily CR Across Locations by Day",
       x = "Day",
       y = "Total Caribou CR") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")

Based on the daily detections it looks like caribou detection rate could be broken down into the following biologically relevant seasons. These dates roughly conform to date ranges for the Beverly, Bathurst, and Ahiak herds listed in various GN and GNWT reports.

Typical caribou seasons (Nagy 2011 and GN use 9 seasons): July 29, 2021 - Oct 20, 2021: Summer/Late Summer Oct 21, 2021 - Oct 31, 2021: Pre-breeding Migration Nov 1, 2021 - Dec 14, 2021: Rut/Post-breeding Migration Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - May 14, 2022: Spring Migration May 15, 2022 - early June, 2022: Calving early June, 2022 - July 30, 2022: Post-calving/Summer July 31, 2022 - Aug 24, 2022: Late Summer

Simplified 4 seasons (Bathurst Range Plan also uses 4): July 29, 2021 - Oct 20, 2021: Summer Oct 21, 2021 - Dec 14, 2021: Fall Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer

##Add simplified 4 seasons

# Define function to assign seasons based on date
assign_season <- function(date) {
  case_when(
    date >= as.Date("2021-07-29") & date <= as.Date("2021-10-20") ~ "Summer",
    date >= as.Date("2021-10-21") & date <= as.Date("2021-12-14") ~ "Fall",
    date >= as.Date("2021-12-15") & date <= as.Date("2022-04-04") ~ "Winter",
    date >= as.Date("2022-04-05") & date <= as.Date("2022-06-21") ~ "Spring",
    date >= as.Date("2022-06-22") & date <= as.Date("2022-08-24") ~ "Summer",
    TRUE ~ "Unknown" # Default case
  )
}

# Create new column "four_season"
mod_dat_daily <- mod_dat_daily %>%
  mutate(four_season = assign_season(day))

# Get the column names
column_names <- names(mod_dat_daily)

# Reorder the columns
new_column_order <- c(column_names[1:4], "snow_season", "four_season", column_names[5:length(column_names)])

# Create the reordered data frame
mod_dat_daily <- mod_dat_daily[, new_column_order]

# Make sure it's a factor
mod_dat_daily <- mod_dat_daily %>% 
            mutate_if(is.character,as.factor)

let’s plot it out

# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
  group_by(year = year(day), day, four_season) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.
# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Locations by Day",
       x = "Day",
       y = "Total Caribou Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")

# Plot the data for capture rate
daily_sum$total_daily_cr <- (daily_sum$total_caribou / daily_sum$total_effort)

ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Daily CR Across Locations by Day",
       x = "Day",
       y = "Total Caribou CR") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")

###Updated colours for thesis

# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
  group_by(year = year(day), day, four_season) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.
# Define the colors for each season using specific colors from Paired
season_colors <- c("Winter" = "#A6CEE3",  # Light Blue
                   "Summer" = "#B2DF8A",  # Light Green
                   "Spring" = "#CAB2D6",  # Light Purple
                   "Fall" = "#FDBF6F")    # Light Orange

# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(fill = "Season", 
       x = "Day",
       y = "Total Caribou Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  scale_fill_manual(values = season_colors)

# Calculate the total daily capture rate
daily_sum$total_daily_cr <- daily_sum$total_caribou / daily_sum$total_effort

# Plot the data for capture rate
ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(fill = "Season",
       x = "Day",
       y = "Caribou detection rate") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  scale_fill_manual(values = season_colors)

##Summary Stats

###Capture Rate by Season

# Summarize the total capture rate by season
seasonal_capture_rate <- daily_sum %>%
  group_by(four_season) %>%
  summarise(total_daily_cr = sum(total_daily_cr))

print(seasonal_capture_rate)
## # A tibble: 4 × 2
##   four_season total_daily_cr
##   <fct>                <dbl>
## 1 Fall                 7.28 
## 2 Spring               0.421
## 3 Summer               1.32 
## 4 Winter               1.77

###Effort by season

# Summarize the total effort by season
seasonal_effort <- daily_sum %>%
  group_by(four_season) %>%
  summarise(total_effort = sum(total_effort))

print(seasonal_effort)
## # A tibble: 4 × 2
##   four_season total_effort
##   <fct>              <int>
## 1 Fall               12514
## 2 Spring             21861
## 3 Summer             30992
## 4 Winter             18460
#calculate mean daily effort by season
seasonal_mean_effort <- mod_dat_daily %>% 
  group_by(four_season) %>% 
  summarise(effort_Mean = mean(n_days_effort)) 

print(seasonal_mean_effort)
## # A tibble: 4 × 2
##   four_season effort_Mean
##   <fct>             <dbl>
## 1 Fall                  1
## 2 Spring                1
## 3 Summer                1
## 4 Winter                1

###Plot daily effort

# Summarize total days effort for each week
total_days_effort <- mod_dat_daily %>%
  group_by(day, four_season) %>%
  summarise(total_effort = sum(n_days_effort), .groups = 'drop')

print(total_days_effort)
## # A tibble: 392 × 3
##    day        four_season total_effort
##    <date>     <fct>              <int>
##  1 2021-07-29 Summer                 6
##  2 2021-07-30 Summer                 6
##  3 2021-07-31 Summer                 6
##  4 2021-08-01 Summer                 6
##  5 2021-08-02 Summer                 6
##  6 2021-08-03 Summer                 6
##  7 2021-08-04 Summer                 7
##  8 2021-08-05 Summer                 7
##  9 2021-08-06 Summer                 7
## 10 2021-08-07 Summer                 7
## # ℹ 382 more rows
# Create the bar plot with colors based on the 'four_season' column
ggplot(total_days_effort, aes(x = day, y = total_effort, fill = four_season)) +
  geom_bar(stat = "identity") +  # Bar plot with black outlines
  scale_fill_manual(values = season_colors) +     # Apply the season colors
  labs(fill = "Season",
       x = "Day of Year",
       y = "Total Effort (Camera Days)") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

#try out a line plot
ggplot(total_days_effort, aes(x = day, y = total_effort)) +
  geom_line(aes(color=four_season)) +  # Line plot
  geom_point(aes(shape=four_season, color=four_season), size=2) +  # Optional: Add points to the line plot
  scale_color_manual(values = season_colors) +
  labs(color = "Season",
       shape = "Season",
       x = "Day",
       y = "Total Effort (Days)") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Multiseason, WEEKLY detections

# Import the total observations dataset
weekly_obs <- read.csv("Processed_Data/wildRtrax/TDN_cam_summaries_week.csv", header=T)

Join this dataframe with the location data, as before:

mod_dat_wk <- left_join(weekly_obs, z_locs)
## Joining with `by = join_by(location)`
# Remove rows with year 2021 and weeks 30, 31, 32, or 33 due to low survey effort
mod_dat_wk <- mod_dat_wk %>%
  filter(!(year == 2021 & week %in% c(30, 31, 32, 33)))

Raw data check:

boxplot(mod_dat_wk$Caribou~mod_dat_wk$dom_habitat_lcc_hexagon,
        las=2,
        xlab="lcc landcover",
        ylab="Habitat use")

Add snow based seasons

# Lets create a new column and give it the value summer
mod_dat_wk$snow_season <- "summer"
mod_dat_wk$snow_season[mod_dat_wk$week %in% c(44:53,1:17)] <- "winter"

# Get the column names
column_names <- names(mod_dat_wk)

# Reorder the columns
new_column_order <- c(column_names[1:3], "snow_season", column_names[4:length(column_names)])

# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]

# make it a factor factors
mod_dat_wk <- mod_dat_wk %>% 
            mutate_if(is.character,as.factor)

Plot caribou detections and rate summed across locations per week

# Convert 'week' and 'year' to a Date object
mod_dat_wk$week_year <- as.Date(paste(mod_dat_wk$year, mod_dat_wk$week, "1", sep="-"), format="%Y-%U-%u")
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid

## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Caribou),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Caribou Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Caribou CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

let’s repeat this process with the new 4-season caribou based dates from the daily detection data

Should the 4-seasons be based now on week instead of the daily detections even though it’s less detailed?

based on daily detection rate: July 29, 2021 - Oct 20, 2021: Summer Oct 21, 2021 - Dec 14, 2021: Fall Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer

based on weekly detection rates? July 26, 2021 - Oct 18, 2021: Summer Oct 25, 2021 - Dec 13, 2021: Fall Dec 20, 2021 - Apr 4, 2022: Winter Apr 11, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer

Add Caribou-centric Seasons

# Let's create a new column for the four seasons based on the specified date ranges
mod_dat_wk$four_season <- "Summer"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-07-29") & mod_dat_wk$week_year <= as.Date("2021-10-20")] <- "Summer"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-10-21") & mod_dat_wk$week_year <= as.Date("2021-12-14")] <- "Fall"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-12-15") & mod_dat_wk$week_year <= as.Date("2022-04-04")] <- "Winter"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2022-04-05") & mod_dat_wk$week_year <= as.Date("2022-06-21")] <- "Spring"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2022-06-22") & mod_dat_wk$week_year <= as.Date("2022-08-24")] <- "Summer"

# Get the column names
column_names <- names(mod_dat_wk)

# Reorder the columns
new_column_order <- c(column_names[1:3], "four_season", column_names[4:length(column_names)])

# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]

# Make it a factor
mod_dat_wk <- mod_dat_wk %>% 
            mutate_if(is.character, as.factor)

Plot it out

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate per 7 days
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Caribou Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Caribou CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

##Summary Stats

###Capture Rate by Season

# Summarize the total capture rate by season
seasonal_capture_rate <- weekly_sum %>%
  group_by(four_season) %>%
  summarise(total_weekly_cr = sum(total_weekly_cr))

print(seasonal_capture_rate)
## # A tibble: 4 × 2
##   four_season total_weekly_cr
##   <fct>                 <dbl>
## 1 Fall                  7.31 
## 2 Spring                0.428
## 3 Summer                1.66 
## 4 Winter                1.69
# Calculate the average weekly capture rate by season
avg_weekly_cr_season <- weekly_sum %>%
  group_by(four_season) %>%
  summarise(avg_weekly_cr = mean(total_weekly_cr))

print(avg_weekly_cr_season)
## # A tibble: 4 × 2
##   four_season avg_weekly_cr
##   <fct>               <dbl>
## 1 Fall               0.914 
## 2 Spring             0.0389
## 3 Summer             0.0875
## 4 Winter             0.106

###Effort by season

# Summarize the total effort by season
seasonal_effort_wk <- weekly_sum %>%
  group_by(four_season) %>%
  summarise(total_effort = sum(total_effort))

print(seasonal_effort_wk)
## # A tibble: 4 × 2
##   four_season total_effort
##   <fct>              <int>
## 1 Fall               12490
## 2 Spring             21579
## 3 Summer             30002
## 4 Winter             19238
#calculate mean weekly effort by season
seasonal_mean_effort_wk <- mod_dat_wk %>% 
  group_by(four_season) %>% 
  summarise(effort_Mean = mean(n_days_effort)) 

print(seasonal_mean_effort_wk)
## # A tibble: 4 × 2
##   four_season effort_Mean
##   <fct>             <dbl>
## 1 Fall               6.21
## 2 Spring             6.84
## 3 Summer             6.49
## 4 Winter             6.21

###Effort by Month

# Let's create a new column for the month/year based on the specified date ranges
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-07-29") & mod_dat_wk$week_year <= as.Date("2021-08-31")] <- "2021-08"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-09-01") & mod_dat_wk$week_year <= as.Date("2021-09-30")] <- "2021-09"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-10-01") & mod_dat_wk$week_year <= as.Date("2021-10-31")] <- "2021-10"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-11-01") & mod_dat_wk$week_year <= as.Date("2021-11-30")] <- "2021-11"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-12-01") & mod_dat_wk$week_year <= as.Date("2021-12-31")] <- "2021-12"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-01-01") & mod_dat_wk$week_year <= as.Date("2022-01-31")] <- "2022-01"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-02-01") & mod_dat_wk$week_year <= as.Date("2022-02-28")] <- "2022-02"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-03-01") & mod_dat_wk$week_year <= as.Date("2022-03-31")] <- "2022-03"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-04-01") & mod_dat_wk$week_year <= as.Date("2022-04-30")] <- "2022-04"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-05-01") & mod_dat_wk$week_year <= as.Date("2022-05-31")] <- "2022-05"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-06-01") & mod_dat_wk$week_year <= as.Date("2022-06-30")] <- "2022-06"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-07-01") & mod_dat_wk$week_year <= as.Date("2022-07-31")] <- "2022-07"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-08-01") & mod_dat_wk$week_year <= as.Date("2022-08-31")] <- "2022-08"

# Get the column names
column_names <- names(mod_dat_wk)

# Reorder the columns
new_column_order <- c(column_names[1:5], "month", column_names[6:length(column_names)])

# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]

# Make it a factor
mod_dat_wk <- mod_dat_wk %>% 
            mutate_if(is.character, as.factor)

#calculate mean weekly effort by month/year
monthly_mean_effort_wk <- mod_dat_wk %>% 
  group_by(month) %>% 
  summarise(effort_Mean = mean(n_days_effort)) 

print(monthly_mean_effort_wk)
## # A tibble: 14 × 2
##    month   effort_Mean
##    <fct>         <dbl>
##  1 2021-08        6.20
##  2 2021-09        6.42
##  3 2021-10        6.83
##  4 2021-11        6.06
##  5 2021-12        5.93
##  6 2022-01        5.57
##  7 2022-02        6.20
##  8 2022-03        6.57
##  9 2022-04        6.87
## 10 2022-05        6.74
## 11 2022-06        7.00
## 12 2022-07        7.00
## 13 2022-08        6.16
## 14 <NA>           1

###Capture Rate by Month

# Calculate the average weekly capture rate by month
avg_weekly_cr_month <- mod_dat_wk %>%
  group_by(month) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort)) %>%
  mutate(avg_weekly_cr = (total_caribou / total_effort) * 7)

print(avg_weekly_cr_month)
## # A tibble: 14 × 4
##    month   total_caribou total_effort avg_weekly_cr
##    <fct>           <int>        <int>         <dbl>
##  1 2021-08             3         1792       0.0117 
##  2 2021-09            21         6074       0.0242 
##  3 2021-10           299         8124       0.258  
##  4 2021-11          1164         8058       1.01   
##  5 2021-12           297         4186       0.497  
##  6 2022-01            53         3742       0.0991 
##  7 2022-02            56         4412       0.0888 
##  8 2022-03            14         7252       0.0135 
##  9 2022-04            60         7853       0.0535 
## 10 2022-05            56         9647       0.0406 
## 11 2022-06            11         8094       0.00951
## 12 2022-07            16         8034       0.0139 
## 13 2022-08           294         5946       0.346  
## 14 <NA>                4           95       0.295

###Plot weekly effort

# Summarize total days effort for each week
total_weeks_effort <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_effort = sum(n_days_effort), .groups = 'drop')

# Convert 'week_year' to Date format if not already done
total_weeks_effort$week_year <- as.Date(total_weeks_effort$week_year, format = "%Y-%m-%d")

# Create the bar plot with colors based on the 'four_season' column
ggplot(total_weeks_effort, aes(x = week_year, y = total_effort, fill = four_season)) +
  geom_bar(stat = "identity") +  # Bar plot with black outlines
  scale_fill_manual(values = season_colors) +     # Apply the season colors
  labs(fill = "Season",
       x = "Week of Year",
       y = "Total Effort (Camera Days)") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +  # Format the x-axis for weekly data
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
## Warning: Removed 1 rows containing missing values (`position_stack()`).

##Consolidate Ungulates and Predators

# Create new columns "Ungulates" and "Predators" by summing the respective columns
mod_dat_wk$Ungulates <- mod_dat_wk$Muskox + mod_dat_wk$Moose
mod_dat_wk$Predators <- mod_dat_wk$Gray.Wolf + mod_dat_wk$Grizzly.Bear

##Plotting Create the same plots for each of the other focal species Muskox

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Muskox),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Muskox Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Muskox CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Muskox), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate  
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Muskox Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Muskox CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Moose

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Moose),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Moose Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Moose CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Moose), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Moose Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Moose CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Gray.Wolf

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Gray.Wolf),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Gray.Wolf), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Grizzly.Bear

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Grizzly.Bear),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Grizzly.Bear), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Ungulates

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Ungulates),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Ungulate Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Ungulate Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Ungulate Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Ungulate CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Ungulates), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Ungulate Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Ungulate Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Ungulate Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Ungulate CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Predators

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Predators),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Predator Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Predator Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Predator Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Predator CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Predators), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Predator Count Summed Across Region by Week",
       x = "Week of the Year",
       y = "Total Predator Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Predator Weekly CR Across Region by Week",
       x = "Week of the Year",
       y = "Total Predator CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

map it out

# Aggregate data by location and week to get total caribou detections
#total_caribou_by_location_week <- mod_dat_wk %>%
 # group_by(location, week) %>%
  #summarise(total_caribou = sum(Caribou))

# Merge total caribou detections with latitude and longitude data
#locations_week <- mod_dat_wk %>%
 # distinct(location, latitude, longitude)

#locations_week <- left_join(locations_week, total_caribou_by_location_week, by = "location") %>%
 # mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))

# Create leaflet maps for each week
#maps_folder <- "maps"
#dir.create(maps_folder, showWarnings = FALSE)

#for (week in unique(mod_dat_wk$week)) {
 # map_data_week <- filter(locations_week, week == !!week)
  #map <- leaflet(map_data_week) %>%
   # addTiles() %>%
    #addCircleMarkers(
     # ~longitude, 
      #~latitude, 
      #radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2), 
      #color = ~if_else(has_detections, "blue", "grey"), 
      #popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
    #)
  
  # Save the map as HTML
  #map_file <- file.path(maps_folder, paste0("map_week_", week, ".html"))
#  saveWidget(map, map_file, selfcontained = TRUE)
#}
# Function to read HTML files and convert them to images
#html_to_image <- function(html_file) {
#  webshot::webshot(html_file, file = tempfile(fileext = ".png"))
#}

# Read each HTML file, convert to images, and store them in a list
#map_images <- lapply(list.files("maps", pattern = "*.html", full.names = TRUE), html_to_image)

# Create the animated GIF
#gif_file <- "caribou_detections.gif"
#image_write(image_join(map_images), gif_file)

# Display the GIF
#browseURL(gif_file)

##Target Analysis

# Reorder the levels of the target_cons column
mod_dat_wk$target_cons <- fct_relevel(mod_dat_wk$target_cons, "None")

#check the levels were reordered correctly. None should be first to act as the control feature. 
levels(mod_dat_wk$target_cons)
## [1] "None"      "Esker"     "Shoreline" "Trail"     "Wetland"
# Group by target_cons and count unique locations within each group
unique_locations_by_target <- mod_dat_wk %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_locations_by_target)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     218
## 2 Esker                      7
## 3 Shoreline                 20
## 4 Trail                     27
## 5 Wetland                   35
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk)
summary(region_wk_aov)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## target_cons     4     50   12.62   9.418 1.34e-07 ***
## Residuals   12882  17260    1.34                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk)
## 
## $target_cons
##                          diff         lwr          upr     p adj
## Esker-None         0.33851816  0.07685521  0.600181105 0.0038263
## Shoreline-None     0.05382748 -0.06082899  0.168483953 0.7031152
## Trail-None        -0.15620948 -0.25213604 -0.060282928 0.0000876
## Wetland-None      -0.04684874 -0.13440328  0.040705800 0.5887736
## Shoreline-Esker   -0.28469068 -0.56653803 -0.002843318 0.0463407
## Trail-Esker       -0.49472764 -0.76948887 -0.219966416 0.0000090
## Wetland-Esker     -0.38536690 -0.65731841 -0.113415387 0.0010530
## Trail-Shoreline   -0.21003696 -0.35206633 -0.068007604 0.0005275
## Wetland-Shoreline -0.10067622 -0.23719081  0.035838367 0.2601690
## Wetland-Trail      0.10936074 -0.01184994  0.230571428 0.0995992

perform the same tests on the new target_hydro (wetland shoreline) and target_linear (trail esker) categories. These targets have been separated because their predicted to have similar effects on caribou occurrences.

# Reorder the levels of the target_hydro column
mod_dat_wk$target_hydro <- fct_relevel(mod_dat_wk$target_hydro, "None")

#check the levels were reordered correctly. None should be first to act as the control feature. 
levels(mod_dat_wk$target_hydro)
## [1] "None"      "Shoreline" "Wetland"
#run anova on hydro targets
region_wk_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk)
summary(region_wk_hydro_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro     2      5   2.716   2.022  0.132
## Residuals    12884  17305   1.343
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk)
## 
## $target_hydro
##                          diff         lwr        upr     p adj
## Shoreline-None     0.06728355 -0.03081301 0.16538011 0.2423553
## Wetland-None      -0.03339267 -0.10800855 0.04122321 0.5459357
## Wetland-Shoreline -0.10067622 -0.21811018 0.01675774 0.1099318

Not significant.

# Reorder the levels of the target_linear column
mod_dat_wk$target_linear <- fct_relevel(mod_dat_wk$target_linear, "None")

#check the levels were reordered correctly. None should be first to act as the control feature. 
levels(mod_dat_wk$target_linear)
## [1] "None"  "Esker" "Trail"
#run anova on linear targets
region_wk_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk)
summary(region_wk_linear_aov)
##                  Df Sum Sq Mean Sq F value  Pr(>F)    
## target_linear     2     45   22.41   16.72 5.6e-08 ***
## Residuals     12884  17265    1.34                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk)
## 
## $target_linear
##                   diff        lwr         upr     p adj
## Esker-None   0.3408211  0.1163509  0.56529136 0.0010882
## Trail-None  -0.1539065 -0.2353337 -0.07247934 0.0000282
## Trail-Esker -0.4947276 -0.7308164 -0.25863884 0.0000027
             Df Sum Sq Mean Sq F value   Pr(>F)    

target_linear 2 45 22.351 16.79 5.23e-08 ***

$target_linear diff lwr upr p adj Esker-None 0.3412555 0.1175085 0.56500256 0.0010243 Trail-None -0.1534721 -0.2346369 -0.07230728 0.0000280

Summer, WEEKLY detection data

# Create a new data frame for summer
mod_dat_wk_summer <- mod_dat_wk[mod_dat_wk$four_season == "Summer", ]

let’s try to map er out

# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_summer %>%
  group_by(location) %>%
  summarise(total_caribou = sum(Caribou))

# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_summer %>%
  distinct(location, latitude, longitude)

locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
  mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))

# Create leaflet map
leaflet(locations) %>%
  addTiles() %>%
  addCircleMarkers(
    ~longitude, 
    ~latitude, 
    radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2), 
    color = ~if_else(has_detections, "blue", "grey"), 
    popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
  )

##Target analysis

#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_smr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_summer)
summary(region_wk_smr_aov)
##               Df Sum Sq Mean Sq F value  Pr(>F)   
## target_cons    4     10  2.5275   3.583 0.00636 **
## Residuals   4619   3258  0.7053                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_summer)
## 
## $target_cons
##                          diff           lwr         upr     p adj
## Esker-None         0.06398224 -2.504331e-01  0.37839758 0.9813072
## Shoreline-None     0.13661331  7.223356e-05  0.27315439 0.0498024
## Trail-None        -0.05386288 -1.714791e-01  0.06375331 0.7220456
## Wetland-None      -0.07473195 -1.819719e-01  0.03250796 0.3164885
## Shoreline-Esker    0.07263107 -2.654815e-01  0.41074361 0.9771589
## Trail-Esker       -0.11784512 -4.487684e-01  0.21307812 0.8679521
## Wetland-Esker     -0.13871419 -4.660932e-01  0.18866483 0.7762729
## Trail-Shoreline   -0.19047619 -3.616388e-01 -0.01931355 0.0203240
## Wetland-Shoreline -0.21134526 -3.755509e-01 -0.04713967 0.0040974
## Wetland-Trail     -0.02086907 -1.697099e-01  0.12797172 0.9954618
#run anova on hydro targets
region_wk_smr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_summer)
summary(region_wk_smr_hydro_aov)
##                Df Sum Sq Mean Sq F value  Pr(>F)   
## target_hydro    2      9   4.373     6.2 0.00205 **
## Residuals    4621   3259   0.705                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_summer)
## 
## $target_hydro
##                          diff         lwr         upr     p adj
## Shoreline-None     0.14180405  0.02513743  0.25847067 0.0122279
## Wetland-None      -0.06954121 -0.16085804  0.02177562 0.1745725
## Wetland-Shoreline -0.21134526 -0.35241929 -0.07027123 0.0013034
            Df Sum Sq Mean Sq F value  Pr(>F)   

target_hydro 2 9 4.376 6.301 0.00185 **

                     diff         lwr         upr     p adj

Shoreline-None 0.14233231 0.02655834 0.25810628 0.0110639

#run anova on linear targets
region_wk_smr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_summer)
summary(region_wk_smr_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2      1  0.7016   0.993  0.371
## Residuals     4621   3267  0.7069
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_summer)
## 
## $target_linear
##                    diff        lwr        upr     p adj
## Esker-None   0.06338601 -0.2066014 0.33337343 0.8462692
## Trail-None  -0.05445910 -0.1544370 0.04551876 0.4082263
## Trail-Esker -0.11784512 -0.4024715 0.16678121 0.5954731

Not significant

Fall, WEEKLY detection data

# Create a new data frame for winter
mod_dat_wk_fall <- mod_dat_wk[mod_dat_wk$four_season == "Fall", ]

plot it out

# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_fall %>%
  group_by(location) %>%
  summarise(total_caribou = sum(Caribou))

# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_fall %>%
  distinct(location, latitude, longitude)

locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
  mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))

# Create leaflet map
leaflet(locations) %>%
  addTiles() %>%
  addCircleMarkers(
    ~longitude, 
    ~latitude, 
    radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2), 
    color = ~if_else(has_detections, "blue", "grey"), 
    popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
  )

##Target analysis

#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_fll_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_fall)
summary(region_wk_fll_aov)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## target_cons    4    126  31.513   5.409 0.000249 ***
## Residuals   2007  11694   5.826                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_fall)
## 
## $target_cons
##                         diff          lwr         upr     p adj
## Esker-None         0.9786811 -0.301492928  2.25885517 0.2259125
## Shoreline-None     0.1107002 -0.506163770  0.72756424 0.9883063
## Trail-None        -0.7448285 -1.257196295 -0.23246068 0.0007122
## Wetland-None      -0.1032963 -0.563718174  0.35712566 0.9731459
## Shoreline-Esker   -0.8679809 -2.267564708  0.53160294 0.4383318
## Trail-Esker       -1.7235096 -3.080283150 -0.36673606 0.0048540
## Wetland-Esker     -1.0819774 -2.419999057  0.25604431 0.1771240
## Trail-Shoreline   -0.8555287 -1.618752187 -0.09230526 0.0189988
## Wetland-Shoreline -0.2139965 -0.943364456  0.51537147 0.9303287
## Wetland-Trail      0.6415322 -0.001878368  1.28494283 0.0511013
#run anova on hydro targets
region_wk_fll_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_fall)
summary(region_wk_fll_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2      4   2.087   0.355  0.701
## Residuals    2009  11815   5.881
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_fall)
## 
## $target_hydro
##                          diff        lwr       upr     p adj
## Shoreline-None     0.17914083 -0.3505139 0.7087955 0.7072086
## Wetland-None      -0.03485566 -0.4285429 0.3588316 0.9765087
## Wetland-Shoreline -0.21399649 -0.8435072 0.4155143 0.7047199

Not significant

#run anova on linear targets
region_wk_fll_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_fall)
summary(region_wk_fll_linear_aov)
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## target_linear    2    122   61.02   10.48 2.96e-05 ***
## Residuals     2009  11698    5.82                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_fall)
## 
## $target_linear
##                   diff        lwr        upr     p adj
## Esker-None   0.9847773 -0.1125523  2.0821070 0.0890824
## Trail-None  -0.7387323 -1.1735898 -0.3038747 0.0002068
## Trail-Esker -1.7235096 -2.8886724 -0.5583468 0.0015442
            Df Sum Sq Mean Sq F value   Pr(>F)    

target_linear 2 122 60.97 10.51 2.87e-05 ***

           diff         lwr         upr     p adj

Trail-None -0.7381761 -1.1721914 -0.3041608 0.0002027

Winter, WEEKLY detection data

# Create a new data frame for winter
mod_dat_wk_winter <- mod_dat_wk[mod_dat_wk$four_season == "Winter", ]

plot it out

# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_winter %>%
  group_by(location) %>%
  summarise(total_caribou = sum(Caribou))

# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_winter %>%
  distinct(location, latitude, longitude)

locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
  mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))

# Create leaflet map
leaflet(locations) %>%
  addTiles() %>%
  addCircleMarkers(
    ~longitude, 
    ~latitude, 
    radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2), 
    color = ~if_else(has_detections, "blue", "grey"), 
    popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
  )

##Target analysis

#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_wtr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_winter)
summary(region_wk_wtr_aov)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## target_cons    4    9.6  2.3953   12.29 6.41e-10 ***
## Residuals   3091  602.2  0.1948                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_winter)
## 
## $target_cons
##                           diff         lwr         upr     p adj
## Esker-None         0.448791822  0.25654490  0.64103875 0.0000000
## Shoreline-None    -0.060501372 -0.15145889  0.03045614 0.3646771
## Trail-None        -0.051592794 -0.12328805  0.02010247 0.2839604
## Wetland-None      -0.004043230 -0.07048943  0.06240297 0.9998297
## Shoreline-Esker   -0.509293194 -0.71877657 -0.29980982 0.0000000
## Trail-Esker       -0.500384615 -0.70225147 -0.29851776 0.0000000
## Wetland-Esker     -0.452835052 -0.65289782 -0.25277228 0.0000000
## Trail-Shoreline    0.008908578 -0.10093047  0.11874762 0.9994679
## Wetland-Shoreline  0.056458142 -0.05002896  0.16294524 0.5971136
## Wetland-Trail      0.047549564 -0.04303986  0.13813898 0.6065337
#run anova on hydro targets
region_wk_wtr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_winter)
summary(region_wk_wtr_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2    0.7  0.3300    1.67  0.188
## Residuals    3093  611.1  0.1976
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_winter)
## 
## $target_hydro
##                           diff         lwr        upr     p adj
## Shoreline-None    -0.060971779 -0.13919861 0.01725505 0.1607324
## Wetland-None      -0.004513637 -0.06136034 0.05233307 0.9810722
## Wetland-Shoreline  0.056458142 -0.03567087 0.14858715 0.3220188

Not significant

#run anova on linear targets
region_wk_wtr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_winter)
summary(region_wk_wtr_linear_aov)
##                 Df Sum Sq Mean Sq F value  Pr(>F)    
## target_linear    2    8.9   4.469   22.93 1.3e-10 ***
## Residuals     3093  602.9   0.195                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_winter)
## 
## $target_linear
##                    diff        lwr         upr     p adj
## Esker-None   0.45359758  0.2887224  0.61847281 0.0000000
## Trail-None  -0.04678703 -0.1075309  0.01395679 0.1677186
## Trail-Esker -0.50038462 -0.6738460 -0.32692325 0.0000000
            Df Sum Sq Mean Sq F value   Pr(>F)    

target_linear 2 8.9 4.470 23.12 1.08e-10 ***

               diff        lwr         upr     p adj

Esker-None 0.45396375 0.2897541 0.61817339 0.0000000

Spring, WEEKLY detection data

# Create a new data frame for winter
mod_dat_wk_spri <- mod_dat_wk[mod_dat_wk$four_season == "Spring", ]

plot it out

# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_spri %>%
  group_by(location) %>%
  summarise(total_caribou = sum(Caribou))

# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_spri %>%
  distinct(location, latitude, longitude)

locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
  mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))

# Create leaflet map
leaflet(locations) %>%
  addTiles() %>%
  addCircleMarkers(
    ~longitude, 
    ~latitude, 
    radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2), 
    color = ~if_else(has_detections, "blue", "grey"), 
    popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
  )

##Target analysis

#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_spr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_spri)
summary(region_wk_spr_aov)
##               Df Sum Sq Mean Sq F value Pr(>F)
## target_cons    4    0.5  0.1283    0.61  0.655
## Residuals   3150  662.0  0.2102
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_spri)
## 
## $target_cons
##                          diff         lwr        upr     p adj
## Esker-None        -0.04172156 -0.28394226 0.20049913 0.9900043
## Shoreline-None     0.02200393 -0.06943901 0.11344687 0.9654090
## Trail-None        -0.03130490 -0.10955760 0.04694780 0.8108308
## Wetland-None      -0.01943744 -0.09048961 0.05161473 0.9454063
## Shoreline-Esker    0.06372549 -0.19251186 0.31996284 0.9610529
## Trail-Esker        0.01041667 -0.24141513 0.26224846 0.9999634
## Wetland-Esker      0.02228412 -0.22740404 0.27197228 0.9992236
## Trail-Shoreline   -0.05330882 -0.16780860 0.06119095 0.7092830
## Wetland-Shoreline -0.04144137 -0.15114608 0.06826335 0.8410457
## Wetland-Trail      0.01186746 -0.08711132 0.11084623 0.9975265
#run anova on hydro targets
region_wk_spr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_spri)
summary(region_wk_spr_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2    0.2  0.1118   0.532  0.588
## Residuals    3152  662.3  0.2101
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_spri)
## 
## $target_hydro
##                          diff         lwr        upr     p adj
## Shoreline-None     0.02591685 -0.05224127 0.10407497 0.7168795
## Wetland-None      -0.01552452 -0.07605275 0.04500371 0.8193488
## Wetland-Shoreline -0.04144137 -0.13568022 0.05279749 0.5572306

Not Significant

#run anova on linear targets
region_wk_spr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_spri)
summary(region_wk_spr_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2    0.3  0.1404   0.668  0.513
## Residuals     3152  662.2  0.2101
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_spri)
## 
## $target_linear
##                    diff         lwr        upr     p adj
## Esker-None  -0.04084507 -0.24866735 0.16697721 0.8895553
## Trail-None  -0.03042840 -0.09689405 0.03603725 0.5306441
## Trail-Esker  0.01041667 -0.20590322 0.22673656 0.9929959

Not Significant

Caribou Sites Only

Subset detection data to only the clusters which captured caribou at a least one of their cameras. Must be subset after joining with covariates to filter by location_site. Must decide on response variable time (i.e. week, biweek) before trying this to cut down on time.

I think weekly is probably the way to go. it’s a good balance of high detail without being so coarse as to become irrelevant for species co-occurrences. It also provides adequate survey replicas when running 4 separate seasonal models.

##UPDATE: 2024-03-21 ##############################################################

Realized that you can’t just subset the data based on caribou detections being greater than 1 because you want to keep ALL cameras and ALL sampling periods at clusters with at least one caribou, whereas the original code was only taking the cameras/sampling periods with caribou and dumping everything else. Need to figure out how to keep only those cameras (i.e. location) which had at least 1 caribou detection at the cluster level (i.e. location_site).

##UPDATE: 2024-03-22 ##############################################################

OKAY I THINK I FIXED IT.

weekly caribou data

# Filter rows where "Caribou" is greater than or equal to 1
caribou_data <- mod_dat_wk %>%
  filter(Caribou >= 1) %>%
  dplyr::select(location, location_site) %>%
  distinct()

# Filter mod_dat_wk for rows with locations that have caribou detections or share the same location_site
mod_dat_wk_cbou <- mod_dat_wk %>%
  filter(location %in% caribou_data$location | location_site %in% caribou_data$location_site)

Better to use biological rationale for choosing the local sites rather than a detection based criteria. Proportion of “TUNDRA” vs “TAIGA”? Compare with locations from caribou based data above

# Create a new dataframe with rows where z.TUNDRA > z.TAIGA
mod_dat_wk_tundra <- mod_dat_wk %>%
  filter(z.TUNDRA > z.TAIGA)

Compare with locations from caribou detection filtered data

# Get unique locations in mod_dat_wk_tundra
unique_locations_tundra <- unique(mod_dat_wk_tundra$location)

# Get unique locations in mod_dat_wk_cbou
unique_locations_cbou <- unique(mod_dat_wk_cbou$location)

# Compare unique locations
common_locations <- intersect(unique_locations_tundra, unique_locations_cbou)
only_in_tundra <- setdiff(unique_locations_tundra, unique_locations_cbou)
only_in_cbou <- setdiff(unique_locations_cbou, unique_locations_tundra)

# Output the results
cat("Common locations:", common_locations, "\n")
## Common locations: BIO-TDN-018-01 BIO-TDN-018-02 BIO-TDN-018-03 BIO-TDN-018-04 BIO-TDN-018-05 BIO-TDN-019-01 BIO-TDN-019-02 BIO-TDN-019-03 BIO-TDN-019-04 BIO-TDN-019-05 BIO-TDN-028-01 BIO-TDN-028-02 BIO-TDN-028-03 BIO-TDN-028-05 BIO-TDN-028-06 BIO-TDN-031-01 BIO-TDN-031-02 BIO-TDN-031-03 BIO-TDN-031-04 BIO-TDN-031-05 BIO-TDN-034-01 BIO-TDN-034-02 BIO-TDN-034-03 BIO-TDN-034-04 BIO-TDN-034-05 BIO-TDN-034-NR BIO-TDN-036-01 BIO-TDN-036-02 BIO-TDN-036-03 BIO-TDN-036-04 BIO-TDN-036-05 BIO-TDN-037-01 BIO-TDN-037-03 BIO-TDN-037-04 BIO-TDN-037-05 BIO-TDN-037-06 BIO-TDN-042-01 BIO-TDN-042-02 BIO-TDN-042-03 BIO-TDN-042-06 BIO-TDN-042-07 BIO-TDN-044-01 BIO-TDN-044-02 BIO-TDN-044-03 BIO-TDN-044-07 BIO-TDN-044-08 BIO-TDN-045-01 BIO-TDN-045-02 BIO-TDN-045-03 BIO-TDN-045-05 BIO-TDN-045-06 BIO-TDN-046-01 BIO-TDN-046-02 BIO-TDN-046-03 BIO-TDN-046-04 BIO-TDN-046-05 BIO-TDN-046-NR BIO-TDN-049-01 BIO-TDN-049-02 BIO-TDN-049-03 BIO-TDN-049-04 BIO-TDN-049-05 BIO-TDN-051-01 BIO-TDN-051-02 BIO-TDN-051-03 BIO-TDN-051-05 BIO-TDN-051-06 BIO-TDN-052-02 BIO-TDN-052-04 BIO-TDN-052-05 BIO-TDN-052-06 BIO-TDN-052-07 BIO-TDN-053-01 BIO-TDN-053-02 BIO-TDN-053-03 BIO-TDN-053-04 BIO-TDN-053-05 BIO-TDN-054-01 BIO-TDN-054-02 BIO-TDN-054-04 BIO-TDN-054-05 BIO-TDN-054-08 BIO-TDN-055-01 BIO-TDN-055-02 BIO-TDN-055-04 BIO-TDN-055-05 BIO-TDN-055-06 BIO-TDN-056-01 BIO-TDN-056-02 BIO-TDN-056-03 BIO-TDN-056-04 BIO-TDN-056-05 BIO-TDN-057-01 BIO-TDN-057-02 BIO-TDN-057-03 BIO-TDN-057-05 BIO-TDN-057-06 BIO-TDN-058-01 BIO-TDN-058-02 BIO-TDN-058-03 BIO-TDN-058-05 BIO-TDN-058-06 BIO-TDN-059-01 BIO-TDN-059-02 BIO-TDN-059-03 BIO-TDN-059-04 BIO-TDN-059-05 BIO-TDN-060-01 BIO-TDN-060-02 BIO-TDN-060-04 BIO-TDN-060-05 BIO-TDN-060-06 BIO-TDN-061-01 BIO-TDN-061-02 BIO-TDN-061-03 BIO-TDN-061-05 BIO-TDN-061-06 BIO-TDN-062-01 BIO-TDN-062-02 BIO-TDN-062-03 BIO-TDN-062-04 BIO-TDN-062-05 BIO-TDN-164-01 BIO-TDN-164-02 BIO-TDN-164-04 BIO-TDN-164-05 BIO-TDN-164-06 BIO-TDN-165-01 BIO-TDN-165-02 BIO-TDN-165-04 BIO-TDN-165-05 BIO-TDN-165-06 BIO-TDN-165-NR BMS-CRU-004-01 BMS-CRU-004-02 BMS-CRU-004-03 BMS-CRU-004-04 BMS-CRU-004-05 BMS-CRU-007-02 BMS-CRU-007-03 BMS-CRU-007-05 BMS-CRU-007-06 BMS-CRU-007-08 BMS-CRU-008-01 BMS-CRU-008-02 BMS-CRU-008-03 BMS-CRU-008-04 BMS-CRU-008-05
cat("Locations only in mod_dat_wk_tundra:", only_in_tundra, "\n")
## Locations only in mod_dat_wk_tundra: BIO-TDN-009-01
cat("Locations only in mod_dat_wk_cbou:", only_in_cbou, "\n")
## Locations only in mod_dat_wk_cbou: BIO-TDN-050-01 BIO-TDN-050-03 BIO-TDN-050-04 BIO-TDN-050-06 BIO-TDN-050-07

Map it out to see the cameras visually

# Create a leaflet map
map <- leaflet() %>%
  addTiles() %>%
  addCircleMarkers(data = mod_dat_wk_tundra[mod_dat_wk_tundra$location %in% common_locations,],
                   lat = ~latitude,
                   lng = ~longitude,
                   label = ~paste(location, " (Common)"),
                   color = "blue",
                   radius = 5) %>%
  addCircleMarkers(data = mod_dat_wk_cbou[mod_dat_wk_cbou$location %in% common_locations,],
                   lat = ~latitude,
                   lng = ~longitude,
                   label = ~paste(location, " (Common)"),
                   color = "red",
                   radius = 5) %>%
  addCircleMarkers(data = mod_dat_wk_tundra[mod_dat_wk_tundra$location %in% unique_locations_tundra,],
                   lat = ~latitude,
                   lng = ~longitude,
                   label = ~paste(location, " (Tundra)"),
                   color = "green",
                   radius = 5) %>%
  addCircleMarkers(data = mod_dat_wk_cbou[mod_dat_wk_cbou$location %in% unique_locations_cbou,],
                   lat = ~latitude,
                   lng = ~longitude,
                   label = ~paste(location, " (Caribou)"),
                   color = "orange",
                   radius = 5) %>%
  addLegend("bottomright",
            colors = c("blue", "red", "green", "orange"),
            labels = c("Common (Tundra)", "Common (Caribou)", "Unique (Tundra)", "Unique (Caribou)"))

# Print the map
map

Very annoyingly, there is one single camera with caribou detections which is not in a majority tundra-type landcover area.

BIO-TDN-050-06 has only two independent detections, both on 2021-11-24, first of a large bachelor group during the day and next of two females that night. Right in the transitional area between tundra and taiga but the camera micro site is very much needle leaf forest.

##Plotting

Local Scale Caribou

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Caribou),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Caribou Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Caribou CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Caribou Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Caribou Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Caribou CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Muskox

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Muskox),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Muskox Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Muskox CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Muskox), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Muskox Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Muskox Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Muskox CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Moose

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Moose),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Moose Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Moose CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Moose), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Moose Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Moose Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Moose CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Gray.Wolf

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Gray.Wolf),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Gray.Wolf), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Gray Wolf Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Gray Wolf CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Grizzly.Bear

# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, snow_season) %>%
  summarise(total_caribou = sum(Grizzly.Bear),total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
#create capture rate 
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)

# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
  group_by(week_year, four_season) %>%
  summarise(total_caribou = sum(Grizzly.Bear), total_effort = sum(n_days_effort))
## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.
# Create capture rate 
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7

# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Count Summed Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear Count") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Total Grizzly Bear Weekly CR Across Local Scale by Week",
       x = "Week of the Year",
       y = "Total Grizzly Bear CR") +
  theme_minimal() +
  scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) 
## Warning: Removed 1 rows containing missing values (`geom_bar()`).

Local, Multiseason, WEEK

##Target analysis

# Reorder the levels of the target_cons column
mod_dat_wk_cbou$target_cons <- fct_relevel(mod_dat_wk_cbou$target_cons, "None")

#check the levels were reordered correctly. None should be first to act as the control feature. 
levels(mod_dat_wk_cbou$target_cons)
## [1] "None"      "Esker"     "Shoreline" "Trail"     "Wetland"
# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target <- mod_dat_wk_cbou %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_cbou_locations_by_target)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     108
## 2 Esker                      7
## 3 Shoreline                 13
## 4 Trail                      5
## 5 Wetland                   20
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)

region_wk_cbou_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou)
summary(region_wk_cbou_aov)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## target_cons    4     26   6.576   2.479  0.042 *
## Residuals   6351  16844   2.652                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou)
## 
## $target_cons
##                          diff        lwr        upr     p adj
## Esker-None         0.14286874 -0.2282522 0.51398970 0.8317097
## Shoreline-None     0.02273957 -0.1869334 0.23241255 0.9983316
## Trail-None        -0.19579310 -0.4866200 0.09503386 0.3522669
## Wetland-None      -0.13667776 -0.2998466 0.02649109 0.1496063
## Shoreline-Esker   -0.12012917 -0.5361447 0.29588642 0.9342055
## Trail-Esker       -0.33866183 -0.8009242 0.12360057 0.2665033
## Wetland-Esker     -0.27954650 -0.6741733 0.11508027 0.2998968
## Trail-Shoreline   -0.21853267 -0.5648252 0.12775985 0.4204316
## Wetland-Shoreline -0.15941733 -0.4083395 0.08950488 0.4048496
## Wetland-Trail      0.05911534 -0.2611655 0.37939619 0.9870349
#run anova on hydro targets
region_wk_cbou_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou)
summary(region_wk_cbou_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)  
## target_hydro    2     14   6.999   2.638 0.0716 .
## Residuals    6353  16857   2.653                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou)
## 
## $target_hydro
##                          diff        lwr         upr     p adj
## Shoreline-None     0.02817302 -0.1513042 0.207650247 0.9280829
## Wetland-None      -0.13124431 -0.2705580 0.008069353 0.0698140
## Wetland-Shoreline -0.15941733 -0.3733242 0.054489516 0.1878788

Not significant (wetland-none is close though at 0.07)

#run anova on linear targets
region_wk_cbou_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou)
summary(region_wk_cbou_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2     12   5.765   2.172  0.114
## Residuals     6353  16859   2.654
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou)
## 
## $target_linear
##                   diff        lwr        upr     p adj
## Esker-None   0.1612578 -0.1565214 0.47903699 0.4594069
## Trail-None  -0.1774040 -0.4258569 0.07104883 0.2153357
## Trail-Esker -0.3386618 -0.7359278 0.05860418 0.1126497

Not Significant

Local, Summer, WEEKLY detection data

# Create a new data frame for summer
mod_dat_wk_cbou_summer <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Summer", ]

##Target analysis

# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_smr <- mod_dat_wk_cbou_summer %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_cbou_locations_by_target_smr)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     108
## 2 Esker                      7
## 3 Shoreline                 13
## 4 Trail                      5
## 5 Wetland                   20
#run anova to see if target is important
region_wk_cbou_smr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_aov)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## target_cons    4     14   3.479   2.453 0.0441 *
## Residuals   2272   3223   1.419                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_smr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_summer)
## 
## $target_cons
##                           diff         lwr         upr     p adj
## Esker-None        -0.021674699 -0.47140773  0.42805833 0.9999326
## Shoreline-None     0.188071890 -0.06114205  0.43728583 0.2379321
## Trail-None        -0.025378402 -0.37742875  0.32667194 0.9996660
## Wetland-None      -0.153482977 -0.35599721  0.04903125 0.2339582
## Shoreline-Esker    0.209746589 -0.29170087  0.71119404 0.7841851
## Trail-Esker       -0.003703704 -0.56341917  0.55601177 1.0000000
## Wetland-Esker     -0.131808279 -0.61176033  0.34814378 0.9446105
## Trail-Shoreline   -0.213450292 -0.62953826  0.20263768 0.6274293
## Wetland-Shoreline -0.341554868 -0.64189096 -0.04121878 0.0165115
## Wetland-Trail     -0.128104575 -0.51801947  0.26181032 0.8981167
#run anova on hydro targets
region_wk_cbou_smr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_hydro_aov)
##                Df Sum Sq Mean Sq F value  Pr(>F)   
## target_hydro    2     14   6.920   4.882 0.00766 **
## Residuals    2274   3223   1.417                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_smr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_summer)
## 
## $target_hydro
##                         diff         lwr         upr     p adj
## Shoreline-None     0.1900115 -0.02308802  0.40311107 0.0918697
## Wetland-None      -0.1515433 -0.32433252  0.02124583 0.0992151
## Wetland-Shoreline -0.3415549 -0.59945484 -0.08365489 0.0054503
           Df Sum Sq Mean Sq F value  Pr(>F)   

target_hydro 2 14 6.912 4.952 0.00715 **

                    diff         lwr         upr     p adj

Wetland-Shoreline -0.3415549 -0.59747378 -0.08563595 0.0050335

******Not significantly different from NONE

#run anova on linear targets
region_wk_cbou_smr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2      0  0.0239   0.017  0.983
## Residuals     2274   3237  1.4236
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_smr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_summer)
## 
## $target_linear
##                     diff        lwr       upr     p adj
## Esker-None  -0.016408814 -0.4019824 0.3691648 0.9945233
## Trail-None  -0.020112518 -0.3212246 0.2809996 0.9865630
## Trail-Esker -0.003703704 -0.4853612 0.4779538 0.9998207

Not Significant

Local, Fall, WEEKLY detection data

# Create a new data frame for winter
mod_dat_wk_cbou_fall <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Fall", ]

##Target analysis

# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_fall <- mod_dat_wk_cbou_fall %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_cbou_locations_by_target_fall)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     104
## 2 Esker                      4
## 3 Shoreline                 13
## 4 Trail                      5
## 5 Wetland                   20
#run anova to see if target is important
region_wk_cbou_fll_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_aov)
##               Df Sum Sq Mean Sq F value Pr(>F)
## target_cons    4     58   14.49   1.403  0.231
## Residuals   1014  10470   10.33
# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_fll_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_fall)
## 
## $target_cons
##                         diff        lwr       upr     p adj
## Esker-None         0.1494062 -1.5712267 1.8700391 0.9992990
## Shoreline-None    -0.1383431 -1.1839603 0.9072742 0.9963502
## Trail-None        -1.0870610 -2.5299283 0.3558062 0.2390066
## Wetland-None      -0.3787046 -1.1907952 0.4333859 0.7071021
## Shoreline-Esker   -0.2877493 -2.2484541 1.6729555 0.9945514
## Trail-Esker       -1.2364672 -3.4348561 0.9619217 0.5385550
## Wetland-Esker     -0.5281108 -2.3748753 1.3186536 0.9359919
## Trail-Shoreline   -0.9487179 -2.6708247 0.7733888 0.5591328
## Wetland-Shoreline -0.2403616 -1.4826485 1.0019253 0.9844157
## Wetland-Trail      0.7083564 -0.8828153 2.2995281 0.7417760
#run anova on hydro targets
region_wk_cbou_fll_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2     13   6.547   0.633  0.531
## Residuals    1016  10515  10.349
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_fll_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_fall)
## 
## $target_hydro
##                         diff        lwr       upr     p adj
## Shoreline-None    -0.0905109 -0.9860893 0.8050675 0.9694576
## Wetland-None      -0.3308725 -1.0246129 0.3628680 0.5022606
## Wetland-Shoreline -0.2403616 -1.3086092 0.8278861 0.8575650

Not Significant

#run anova on linear targets
region_wk_cbou_fll_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2     41   20.37   1.974  0.139
## Residuals     1016  10487   10.32
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_fll_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_fall)
## 
## $target_linear
##                   diff       lwr       upr     p adj
## Esker-None   0.2159652 -1.255703 1.6876329 0.9367034
## Trail-None  -1.0205021 -2.252478 0.2114736 0.1269366
## Trail-Esker -1.2364672 -3.124384 0.6514495 0.2739310

Not Significant

Local, Winter, WEEKLY detection data

# Create a new data frame for winter
mod_dat_wk_cbou_winter <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Winter", ]

##Target analysis

# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_wtr <- mod_dat_wk_cbou_winter %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_cbou_locations_by_target_wtr)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     101
## 2 Esker                      3
## 3 Shoreline                 13
## 4 Trail                      5
## 5 Wetland                   19
#run anova to see if target is important
region_wk_cbou_wtr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_aov)
##               Df Sum Sq Mean Sq F value  Pr(>F)    
## target_cons    4    7.3  1.8336   4.833 0.00071 ***
## Residuals   1551  588.4  0.3794                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_wtr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_winter)
## 
## $target_cons
##                           diff        lwr         upr     p adj
## Esker-None         0.378046595  0.1073555  0.64873769 0.0013369
## Shoreline-None    -0.117827191 -0.2910513  0.05539690 0.3409280
## Trail-None        -0.017921147 -0.2374046  0.20156234 0.9994520
## Wetland-None      -0.027804469 -0.1485355  0.09292653 0.9704193
## Shoreline-Esker   -0.495873786 -0.8092576 -0.18248998 0.0001603
## Trail-Esker       -0.395967742 -0.7371068 -0.05482873 0.0134634
## Wetland-Esker     -0.405851064 -0.6935641 -0.11813801 0.0011496
## Trail-Shoreline    0.099906044 -0.1704802  0.37029233 0.8512838
## Wetland-Shoreline  0.090022723 -0.1087528  0.28879827 0.7297651
## Wetland-Trail     -0.009883322 -0.2500460  0.23027940 0.9999640
#run anova on hydro targets
region_wk_cbou_wtr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2    1.8  0.8794   2.299  0.101
## Residuals    1553  594.0  0.3825
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_wtr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_winter)
## 
## $target_hydro
##                          diff         lwr        upr     p adj
## Shoreline-None    -0.12933027 -0.27820950 0.01954896 0.1035810
## Wetland-None      -0.03930755 -0.14267890 0.06406380 0.6454035
## Wetland-Shoreline  0.09002272 -0.08142462 0.26147007 0.4345380

Not Significant

#run anova on linear targets
region_wk_cbou_wtr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_linear_aov)
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## target_linear    2    6.0  2.9812    7.85 0.000405 ***
## Residuals     1553  589.8  0.3798                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_wtr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_winter)
## 
## $target_linear
##                     diff        lwr        upr     p adj
## Esker-None   0.390887208  0.1591766  0.6225978 0.0002334
## Trail-None  -0.005080534 -0.1925603  0.1823992 0.9977743
## Trail-Esker -0.395967742 -0.6891631 -0.1027724 0.0044555
            Df Sum Sq Mean Sq F value   Pr(>F)    

target_linear 2 6.0 2.9911 7.94 0.000371 ***

                diff        lwr        upr     p adj

Esker-None 0.391574966 0.1607959 0.6223541 0.0002123

Local, Spring, WEEKLY detection data

# Create a new data frame for spring
mod_dat_wk_cbou_spri <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Spring", ]

##Target analysis

# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_spr <- mod_dat_wk_cbou_spri %>%
  group_by(target_cons) %>%
  summarise(unique_locations = n_distinct(location))

# Print the result
print(unique_cbou_locations_by_target_spr)
## # A tibble: 5 × 2
##   target_cons unique_locations
##   <fct>                  <int>
## 1 None                     102
## 2 Esker                      3
## 3 Shoreline                 13
## 4 Trail                      5
## 5 Wetland                   19
#run anova to see if target is important
region_wk_cbou_spr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_aov)
##               Df Sum Sq Mean Sq F value Pr(>F)
## target_cons    4    0.7  0.1677   0.383  0.821
## Residuals   1499  656.9  0.4382
# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_spr_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_spri)
## 
## $target_cons
##                          diff        lwr        upr     p adj
## Esker-None        -0.08715596 -0.4393827 0.26507074 0.9616486
## Shoreline-None     0.01520624 -0.1543138 0.18472626 0.9992054
## Trail-None        -0.03261051 -0.2824718 0.21725078 0.9965482
## Wetland-None      -0.04813157 -0.1857687 0.08950553 0.8750307
## Shoreline-Esker    0.10236220 -0.2807866 0.48551098 0.9496320
## Trail-Esker        0.05454545 -0.3703035 0.47939438 0.9967615
## Wetland-Esker      0.03902439 -0.3311241 0.40917290 0.9984998
## Trail-Shoreline   -0.04781675 -0.3396557 0.24402224 0.9917116
## Wetland-Shoreline -0.06333781 -0.2675027 0.14082708 0.9156957
## Wetland-Trail     -0.01552106 -0.2900697 0.25902753 0.9998725
#run anova on hydro targets
region_wk_cbou_spr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_hydro_aov)
##                Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro    2    0.4  0.2109   0.482  0.618
## Residuals    1501  657.2  0.4378
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_spr_hydro_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_spri)
## 
## $target_hydro
##                          diff        lwr        upr     p adj
## Shoreline-None     0.01874446 -0.1262726 0.16376147 0.9505772
## Wetland-None      -0.04459336 -0.1621119 0.07292516 0.6465706
## Wetland-Shoreline -0.06333781 -0.2386331 0.11195750 0.6733747

Not Significant

#run anova on linear targets
region_wk_cbou_spr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_linear_aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## target_linear    2    0.2  0.1054   0.241  0.786
## Residuals     1501  657.4  0.4380
# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_spr_linear_aov)

# View the results
print(tukey_results)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_spri)
## 
## $target_linear
##                    diff        lwr       upr     p adj
## Esker-None  -0.08157525 -0.3831898 0.2200393 0.8011177
## Trail-None  -0.02702979 -0.2403878 0.1863282 0.9524751
## Trail-Esker  0.05454545 -0.3102870 0.4193779 0.9344405

Not Significant

#SUMMARIZE

Create tables summarizing effort and independent detections for each focal species in each subset dataset.

# List of dataset names
dataset_names <- c("mod_dat_wk", "mod_dat_wk_summer", "mod_dat_wk_fall", "mod_dat_wk_winter", "mod_dat_wk_spri",
                    "mod_dat_wk_cbou", "mod_dat_wk_cbou_summer", "mod_dat_wk_cbou_fall", "mod_dat_wk_cbou_winter",
                    "mod_dat_wk_cbou_spri")

# Get the dataframes from the environment
datasets <- mget(dataset_names)

# Function to create summary dataframe for each dataset
summarize_dataset <- function(df, dataset_name) {
  summary_df <- data.frame(
    Dataset = dataset_name,
    Total_Effort = sum(df$n_days_effort),
    Caribou_Detections = sum(df$Caribou, na.rm = TRUE),
    Muskox_Detections = sum(df$Muskox, na.rm = TRUE),
    Moose_Detections = sum(df$Moose, na.rm = TRUE),
    Gray_Wolf_Detections = sum(df$Gray.Wolf, na.rm = TRUE),
    Grizzly_Bear_Detections = sum(df$Grizzly.Bear, na.rm = TRUE)
  )
  return(summary_df)
}

# Create summary dataframes for each dataset
summary_dfs <- mapply(summarize_dataset, datasets, dataset_names, SIMPLIFY = FALSE)

# Combine summary dataframes into a single dataframe
summary_all <- do.call(rbind, summary_dfs)

# Print summary dataframe
print(summary_all)
##                                       Dataset Total_Effort Caribou_Detections
## mod_dat_wk                         mod_dat_wk        83309               2348
## mod_dat_wk_summer           mod_dat_wk_summer        30002                372
## mod_dat_wk_fall               mod_dat_wk_fall        12490               1633
## mod_dat_wk_winter           mod_dat_wk_winter        19238                224
## mod_dat_wk_spri               mod_dat_wk_spri        21579                119
## mod_dat_wk_cbou               mod_dat_wk_cbou        40039               2348
## mod_dat_wk_cbou_summer mod_dat_wk_cbou_summer        14733                372
## mod_dat_wk_cbou_fall     mod_dat_wk_cbou_fall         5971               1633
## mod_dat_wk_cbou_winter mod_dat_wk_cbou_winter         9226                224
## mod_dat_wk_cbou_spri     mod_dat_wk_cbou_spri        10109                119
##                        Muskox_Detections Moose_Detections Gray_Wolf_Detections
## mod_dat_wk                           746              180                  130
## mod_dat_wk_summer                    599               97                   42
## mod_dat_wk_fall                       54               13                   34
## mod_dat_wk_winter                     22               22                   15
## mod_dat_wk_spri                       71               48                   39
## mod_dat_wk_cbou                      212               29                   57
## mod_dat_wk_cbou_summer               180               19                    6
## mod_dat_wk_cbou_fall                  11                2                   27
## mod_dat_wk_cbou_winter                 3                1                   13
## mod_dat_wk_cbou_spri                  18                7                   11
##                        Grizzly_Bear_Detections
## mod_dat_wk                                  38
## mod_dat_wk_summer                           25
## mod_dat_wk_fall                              0
## mod_dat_wk_winter                            0
## mod_dat_wk_spri                             13
## mod_dat_wk_cbou                             37
## mod_dat_wk_cbou_summer                      25
## mod_dat_wk_cbou_fall                         0
## mod_dat_wk_cbou_winter                       0
## mod_dat_wk_cbou_spri                        12
#save table
write.csv(summary_all, file = "Tables/GLMM_DataSet_CountandEffort_Comparison.csv", row.names = F)
                                  Dataset Total_Effort Caribou_Detections Muskox_Detections Moose_Detections

mod_dat_wk mod_dat_wk 83309 2348 746 180 mod_dat_wk_summer mod_dat_wk_summer 30002 372 599 97 mod_dat_wk_fall mod_dat_wk_fall 12490 1633 54 13 mod_dat_wk_winter mod_dat_wk_winter 19238 224 22 22 mod_dat_wk_spri mod_dat_wk_spri 21579 119 71 48 mod_dat_wk_cbou mod_dat_wk_cbou 40039 2348 212 29 mod_dat_wk_cbou_summer mod_dat_wk_cbou_summer 14733 372 180 19 mod_dat_wk_cbou_fall mod_dat_wk_cbou_fall 5971 1633 11 2 mod_dat_wk_cbou_winter mod_dat_wk_cbou_winter 9226 224 3 1 mod_dat_wk_cbou_spri mod_dat_wk_cbou_spri 10109 119 18 7 Gray_Wolf_Detections Grizzly_Bear_Detections mod_dat_wk 130 38 mod_dat_wk_summer 42 25 mod_dat_wk_fall 34 0 mod_dat_wk_winter 15 0 mod_dat_wk_spri 39 13 mod_dat_wk_cbou 57 37 mod_dat_wk_cbou_summer 6 25 mod_dat_wk_cbou_fall 27 0 mod_dat_wk_cbou_winter 13 0 mod_dat_wk_cbou_spri 11 12

#SAVE DATA AS CSV

##Regional data sets

#multiseason
#write.csv(mod_dat_wk, file="Processed_Data/glmm_wk_region_multi.csv", row.names = F)
#summer
#write.csv(mod_dat_wk_summer, file="Processed_Data/glmm_wk_region_smr.csv", row.names = F)
#fall
#write.csv(mod_dat_wk_fall, file="Processed_Data/glmm_wk_region_fall.csv", row.names = F)
#winter
#write.csv(mod_dat_wk_winter, file="Processed_Data/glmm_wk_region_wtr.csv", row.names = F)
#spring
#write.csv(mod_dat_wk_spri, file="Processed_Data/glmm_wk_region_spr.csv", row.names = F)

##Local data sets

#multiseason
#write.csv(mod_dat_wk_cbou, file="Processed_Data/glmm_wk_local_multi.csv", row.names = F)
#summer
#write.csv(mod_dat_wk_cbou_summer, file="Processed_Data/glmm_wk_local_smr.csv", row.names = F)
#fall
#write.csv(mod_dat_wk_cbou_fall, file="Processed_Data/glmm_wk_local_fall.csv", row.names = F)
#winter
#write.csv(mod_dat_wk_cbou_winter, file="Processed_Data/glmm_wk_local_wtr.csv", row.names = F)
#spring
#write.csv(mod_dat_wk_cbou_spri, file="Processed_Data/glmm_wk_local_spr.csv", row.names = F)

#regional only dataset

#mod_dat_wk_region_only <- subset(mod_dat_wk, spatial_scale == "regional")